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ABSTRACT 


This  paper  gives  a  general  method  for  the  stable  evaluation  of  inner  products  of 
multivariate  simplex  splines.  The  method  is  based  on  a  recurrence  relation  for  these  inner 


products.  The  base  cases  for  this  recurrence  relation  are  handled  by  triangulating  certain 
simploids  into  simplices.  ■  -s  r  ^ ^  j- '•  7 
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SIGNIFICANCE  AND  EXPLANATION 


For  many  types  of  projection  methods  (Rayleigh-Ritz-Galerkin.  for  example),  it  is 
necessary  to  compute  inner  products  of  the  basis  elements  in  order  to  be  able  to  implement 
the  scheme.  This  paper  discusses  a  method  for  evaluating  these  inner  products  in  the  case 
where  the  basis  elements  are  simplex  splines. 
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The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive  summary  lies 
with  MRC,  and  not  with  the  author  of  this  report. 


THE  EVALUATION  OF  INNER  PRODUCTS  OF 
MULTIVARIATE  SIMPLEX  SPLINES 
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Thomas  A.  Grandine*'^ 

In  many  practical  problems,  it  is  necessary  to  be  able  to  compute  the  value  of  the 
inner  product  of  any  two  simplex  splines.  This  is  typically  the  case  for  many  types  of 
projection  methods,  including  Rayleigh-Ritz-Galerkin  and  least  squares.  While  various 
types  of  quadrature  methods  can  be  used  for  this  purpose,  Dahmen  and  Micchelli  [DM81  j 
report  that  far  better  results  are  obtained  if  exact  values  of  these  inner  products  are  used. 

Developing  an  algorithm  for  evaluating  these  inner  products  is  the  goal  of  this  paper. 

The  multivariate  polyhedral  B-spline  Mb  is  defined  as  a  distribution  by 

[  MBf  :=  f  foP,  (1) 

JR’"  Jb 

where  /  is  an  arbitrary  m-variate  function,  B  is  a  polyhedral  body  in  R",  n  >  m,  and 
P  :R"  R"*  is  a  linear  map.  Typically,  P  is  chosen  to  be  the  canonical  projector,  i.e. 

P  :  X  >-*  x(t)'^,.  This  is  a  generalization  of  an  earlier  definition  due  to  Micchelli  in  which 
B  is  chosen  specifically  to  be  a  simplex.  In  particular,  if  lo,  xi,  •••,  Xn  are  points  in  R”. 
and  if  \A]  denotes  the  convex  hull  of  A,  then  the  multivariate  simplex  spline  M{-  xo. 
is  defined  to  be  [M80j: 


_ 1 

vol„jio. 


(2) 
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This  definition,  which  is  a  convenient  reformulation  of  an  earlier  definition  of  de  Boor 
:B76  ,  is  a  special  case  of  definition  (l)  except  for  the  difference  in  normalization. 

The  purpose  of  this  paper  is  to  describe  a  method  for  evaluating  inner  products  of 
simplex  splines  (with  the  normalization  given  in  (2)).  The  inner  product  referred  to  here 
is  the  one  given  by  the  following  integral: 


(/ 


,9)  :=  / 

Jn." 


fg- 


An  identity  can  be  used  to  evaluate  the  integral  in  the  C2ise  where  /  and  g  happen  to  be 
two  polyhedral  splines,  say  Mb  and  Me,  defined  as  in  (1)  by 

/  Mb/:=  /  foP 

JR"-  JB 

I  Mcf~[foQ, 

Jr."'  Jc 

for  two  ordinary  linear  maps  P  and  Q  into  R*”.  From  [DM82-2],  it  is  known  that  the 
convolution  of  two  polyhedral  splines  satisfies 


Mb^c 


MB(y  -  x)Mc[y)dy, 


(3) 


where  B  ^  C  is  the  polyhedral  set  which  is  the  usual  Cartesian  product  of  the  polyhedral 
sets  B  and  C,  i.e.  B  y  C  {(r,s)|r  G  B,s  G  C}.  Setting  i  =  0  in  this  equation  results 
in  the  identity 


MB(i)Mc(x)rfi  =  Mbxc(O), 


(4) 


where 


MB^c(x)f{T)dx 


Ur 


f{Qy  -  Px)dydx. 


(5) 


If  B  and  C  are  simplices.  then  this  relation  says  that  the  value  of  the  inner  product  of  two 
simplex  splines  is  given  by  the  value  of  a  certain  simploid  spline  at  the  origin.  A  simploid 


is  defined  as  the  Cartesian  product  of  two  simplices.  Thus,  if  an  accurate  way  can  be 
found  to  evaluate  this  simploid  spline,  then  an  accurate  way  has  been  found  to  compute 
the  value  of  the  desired  inner  product. 

Dahmen  and  Micchelli  iDM81-2;  have  proposed  the  use  of  the  recurrence  relation  for 
polyhedral  splines  for  this  purpose.  Specialized  to  the  case  of  the  simploid,  and  then 
rewritten  in  terms  of  inner  products,  this  recurrence  relation  takes  on  the  following  form 
!DM81-2l: 

Theorem  1:  If  B  =  [io,...,in],  C  -  [yo,...,ypli  and  Er=o®'«®»  = 

I  M{x\xo,....,Xn)M{x\yo,....,yp)dx  =  - - - ( 

7r...  n  +  p  -  m  V 

"  /■ 

nVo!,  /  M(i|io,...,Xi_i,ii+i,..Marn)Af(x|yo,...,yp)rfi+ 
i^o  ^  ’ 

f  A/(i(xo,...,x„)Af(x(yo,...,y>-i,y;+i,...,yp)rfx). 

7^0  / 

This  recurrence  relation  can  be  implemented  using  the  linear  programming  approatch 
of  [G84].  This  technique  requires,  at  the  first  step,  a  solution  to  the  following  problem: 


Find  Ot  and  I3j  so  that 


£;o,x,  -  £/?jyy  =0 

|  =  CP  >=0 

±a.  =  l 

t=0 


a,  >0,  i  =  0, ...,  n 


/?!  >0,  y  =  o,...,p. 


As  described  in  G84;,  problems  of  this  sort  ran  be  via  the  dual  simplex  method  for 

linear  programming  problems.  This  method  always  results  in  solutions  for  which  all  but 
m  2  of  the  o,  and  are  zero.  This  is  because  the  problem  (7)  has  only  m  •+  2  equality 
constraints  in  it.  and  solutions  to  linear  programming  problems  can  always  be  made  to 
satisfy  a  complementarity  condition  (see  |Da63'  and  iMaTSj).  This  number  is  a  drastic 
improvement  over  the  2m  -r  2  given  by  Oahmen  and  Micchelli  in  [DMSl]. 

In  any  case,  the  implementation  of  the  recurrence  relation  (6)  proceeds  just  as  in  [G84]. 
The  snag  appears  when  one  or  the  other  (or  both)  of  the  simplex  splines  appearing  on  the 
left-hand-side  of  (6)  are  piecewise  constant  functions,  which  happens  when  a  simplex  spline 
has  exactly  m  +  1  knots  present.  In  this  case,  the  recurrence  relation  cannot  be  applied, 
for  then  some  of  the  functions  appearing  on  the  right-hand-side  of  the  recurrence  relation 
will  not  be  well-defined.  Dahmen  and  Micchelli  conquer  this  obstacle  by  expressing  each 
of  the  simplex  splines  in  the  integral  on  the  left-hand-side  of  (6)  as  a  linear  combination 
of  cone  splines.  Then,  since  the  convolution  of  two  cone  splines  is  again  a  cone  spline, 
they  are  able  to  write  this  inner  product  as  a  linear  combination  of  cone  splines  (see  [D79], 
[D80),  and  (DM81-2]),  each  of  which  may  be  evaluated  using  the  recurrence  relation  for 
cone  splines. 

There  is  a  way,  however,  to  avoid  going  to  the  trouble  of  implementing  the  recurrence 
relation  for  cone  splines.  In  fact,  the  cone  splines  can  be  dispensed  with  altogether  by  not¬ 
ing,  as  was  done  already  in  [H82]  and  |DM82],  that  any  simploid  can  easily  be  triangulated 
merely  by  arranging  its  vertices  in  a  rectangular  grid  and  tracing  the  various  paths  through 
this  grid.  This  makes  it  possible  to  express  any  simploid  spline  as  a  linear  combination  of 
simplex  splines,  each  of  which  may  be  evaluated  using  the  method  described  in  |G84). 


In  what  follows,  let  Ei  :=  [io,...,x„]  and  E2  lyoi-Mj/pl-  With  t  :=  let 

be  a  collection  of  simplices  such  that  {<7,}  is  the  triangulation  of  £]  x  £2  that  is 
produced  via  the  following  construction  iH82i;  Consider  the  following  grid,  representing 
the  set  El  X  E2, 

(xo.Vp)  (xi,yp)  ...  (in,yr) 

(xo,yi)  (xi,yt)  ...  (x„,y,) 

(xo,yo)  (ii,yo)  •••  (xn,yo) 

Take  all  the  paths  through  this  grid  given  by  o  =  {so, •■.•■Sn-4-;,},  where  so  =  (a:o,yo), 
s„+p  =  (x„,yp),  and  if  s,  =  (xj,y<),  then  must  be  either  (xj-uye)  or  (xj,ye^i). 
There  are  (^)  different  paths  through  the  grid.  It  can  be  shown  that  the  set  of  all  paths 
through  the  grid  is  a  triangulation  of  the  simploid  (Lemma  3  in  IH82]). 

This  construction  ensures  the  truth  of 


Lemma  1: 


voloj  = 


VolEi  ^2 


Proof:  Let  Pj  be  the  canonical  projector  from  Ei  x  E2  into  Ej.  Then 

volEi  X  E2  volEi  •  V0IE2 

t  “  i 

_  nlvoIFiOi  •  p!volP20t 
(n  -f  p)! 

=  voloi. 

With  this  lemma,  it  is  possible  to  prove 


Theorem  2: 


f  M(x|Xf„...,x„)Af(x|yo,...,yp)dx  =  (1/t)  V  — . 


Proof:  By  definition, 


so  that  Afj:,  .  Then 


As  it  turns  out,  the  terms  on  the  right-hand-side  of  this  formula  are  all  normalized 
in  the  standard  way  for  simplex  splines.  Hence,  the  simplex  splines  appearing  on  the 
right-hand-side  of  (8)  are  those  obtained  by  considering  all  possible  paths  through  the 
grid; 

Xo-  Vp  Xi-  yp  ...  Xn-  yp 

Xo-  Vl  Xi  -yj  ...  Xn  -  yi 

Xo  ~  Vo  Xi  —  yo  ...  Xn  —  J/O 

Formula  (8)  provides  an  easy  way  of  evaluating  inner  products  when  the  recurrence 
relation  (6)  cannot  be  used  because  one  or  more  of  the  simplex  splines  in  the  integral 
have  only  m  -  1  knots.  It  also  provides  an  easy  way  to  dispense  with  (6)  altogether,  since 
(8)  is  true  independent  of  the  number  of  knots  appearing  in  the  simplex  splines  in  the 


integral.  As  elegant  and  clean  as  this  approach  is,  however,  it  is  much  less  efficient  than 
using  recurrence  relation  (6)  whenever  possible,  and  then  using  relation  (8)  the  rest  of  the 
time. 

This  is  best  seen  by  considering  the  following  calculation.  Suppose  that  the  value  of 
an  inner  product  of  a  simplex  spline  with  n  +  1  knots  and  a  simplex  spline  of  p  +  1  knots 
is  desired.  To  obtain  this  value  using  formula  (8)  will  require  the  evaluation  of 
simplex  splines  of  degree  n  -t-  p  —  m.  Each  of  these  simplex  splines  requires  m  1  times  as 
much  work  to  evaluate  as  does  a  simplex  spline  of  degree  n  +  p  —  m  -  1 .  Thus,  the  total 
amount  of  work  required  is  roughly  equivalent  to  that  of  evaluating  (m  +  1)  simplex 

splines  of  degree  n  +  p  —  m  —  1.  Now  consider  what  happens  when  one  application  of  the 
recurrence  relation  (6)  is  applied,  followed  by  this  triangulation  technique.  The  rcK'urrence 
relation  will  produce  a  sum  of  inner  products,  say  inner  products  of  n  +  1  knot  and 
p  knot  simplex  splines,  and  ^2  inner  products  of  n  knot  and  p  +  J  knot  simplex  splines. 
Here  Si  +  6^  =  m  ■+■  2,  as  remarked  earlier.  After  triangulating  each  of  these  simploids, 
the  total  number  of  degree  n  +  p  —  m  -  1  simplex  splines  which  need  to  be  evaluated  is 
+  ^2  Since  Si  <  m  +  1  and  S2  <  m  +  1,  with  equality  in  at  most 

one  of  them,  it  follows  that  ^ip  +  ^2”  <  (m  +  l)(n  p).  Multiplying  this  inequality  by 
(n  +  p  -  l)!/(n!p!)  yields 

Thus,  it  is  always  more  efficient  to  employ  recurrence  relation  (6)  whenever  that  is  possible. 

Example:  Suppose  that  n  =  p  =  5  and  m  =  2  above.  Then  using  the  recurrence 
relation  (6)  once  requires  evaluation  of  504  simplex  splines  of  degree  7,  while  not  using  it 


at  all  requires  evaluation  of  756  simplex  spiines  of  dp’r«:«  7.  naif  agaiti  as  much  work.  Of 
course,  it  would  make  sense  to  employ  (6)  as  often  as  possible  to  compute  inner  products. 
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